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Abstract 

We consider a queueing system consisting of two non-identical exponential servers, where 
each server has its own dedicated queue and serves the customers in that queue FCFS. 
Customers arrive according to a Poisson process and join the queue promising the shortest 
expected delay, which is a natural and near-optimal policy for systems with non-identical 
servers. This system can be modeled as an inhomogeneous random walk in the quadrant. 

By stretching the boundaries of the compensation approach we prove that the equilibrium 
distribution of this random walk can be expressed as a series of product-forms that can be 
determined recursively. The resulting series expression is directly amenable for numerical 
calculations and it also provides insight in the asymptotic behavior of the equilibrium 
probabilities as one of the state coordinates tends to infinity. 

1 Introduction 

In this paper we analyze the performance of a system with two servers under the shortest 
expected delay (SED) routing policy. This routing policy assigns an arriving customer to the 
queue that has the shortest expected delay (sojourn time), where delay refers to the waiting 
time plus the service time. This policy arises naturally in various application areas, and poses 
considerable mathematical challenges. 

In particular, we focus on a queueing system with two servers, where each server has its 
own queue with unlimited buffer capacity. All service times are independent and exponentially 
distributed, but the two servers have different service rates, i.e. respectively 1 and s. In both 
queues customers are served in FCFS-order. Customers arrive according to a Poisson process 
with rate A and upon arrival join one of the two queues according to the following mechanism; 
Let qi and q 2 be the number of customers in queue 1 and 2, respectively, including a possible 
customer in service. For an arriving customer, the expected delay in queue 1 is (71 -|- 1 and 
in queue 2 is ((72 + 1)/'S- Th® SED routing policy assigns an arriving customer to queue 1 if 
qi + 1 < {q 2 + 1)/s and to queue 2 if (71 -|- 1 > ((72 -t- l)/s. In case the expected delays in both 
queues are equal, i.e. qi + 1 = {q 2 + l)/^, the arriving customer joins queue 1 with probability 
q and queue 2 with probability 1 — q. Once a customer has joined a queue, no switching or 
jockeying is allowed. The service times are assumed independent of the arrival process and the 
customer decisions. This system is stable if and only if, see e.g. m Theorem 1], 

P •= (1 + ■s) < 1- (1-1) 
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^ qi + l = {q 2 + l)/s 


<12 



Figure 1: Transition rate diagram of the Markov process on the state space (gi, 52 ) with s = |. 
Only transitions corresponding to arrivals that reach or cross the dashed line are shown. 


We refer to this specific queueing system as the SED system. The SED system can be modeled 
as a two-dimensional Markov process on the states {qi, ^ 2 ) with the transition rate diagram as 
in Figure From the transition rate diagram it is evident that this state description leads to 
an inhomogeneous random walk in the quadrant, making an exact analysis extremely difficult. 
The inhomogeneous behavior occurs along the line s{qi -|- 1) = ^2 + 1 and it can be expected 
that the solution structure of the stationary probabilities above and below this line will be 
different. Moreover, for 57 ^!, this line divides the state space in unequal proportions, further 
increasing the complexity of an exact analysis. 

Under the assumption of identical service rates, SED routing becomes join the shortest 
queue (JSQ) routing which is known to minimize the mean expected delay [29]. If the service 
rates of the two servers are different, then JSQ is not optimal. As Whitt [29| points out, if 
the system if fully observable and we are a priori knowledgeable of the service times of waiting 
customers, for example, by looking at the shopping baskets in a supermarket, then the natural 
choice is to join the queue promising the smallest sojourn time in the system, instead of the 
shorter of the queues. However, this type of information is too detailed and might not always 
be available. In particular, there are many situations in which we have limited knowledge of the 
service times of the waiting customers. Such systems include the teller waiting lines, production 
facilities, communication networks, etc. If the only available information is that of the number 
of waiting customers at each queue, then it is quite natural, although not necessarily optimal, 
to choose the shortest queue routing |ll[5l[9llini[I3ll8l|2ll[2ai26l[271|30l|3l|. However, if on 
top of the number of waiting customers we can estimate the expected service times at the two 
queues, then the natural choice is to route the customers according to SED, rather than JSQ. 

SED routing seems to be a natural choice, but in practice this choice is not always optimal, 
since it does not minimize the mean stationary delay naEei. Eor the two non-identical server 
setting, Hajek |23j solves a Markov decision problem and proves that the optimal routing 
policy is of the threshold-type. However, SED routing exhibits relatively good performance at 
both ends of the range of system utilizations, but performs slightly worse than other policies 
at medium utilizations, see [7|. Furthermore, Foschini [T9| has shown that SED routing is 
asymptotically optimal and results in complete resource pooling in the heavy traffic limit. 

State-dependent routing policies such as JSQ and SED are typically difficult to analyze. 
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For the stationary behavior of queueing systems with SED routing very little is known. The 
difficulty in analyzing this type of models is evident from the analysis of the JSQ policy for 
two identical parallel exponential servers [22]. The first major steps towards its analysis were 
made in [13 [26], using a uniformization approach that established that the generating function 
of the equilibrium distribution is meromorphic. Thus, a partial fraction expansion of the 
generating function of the joint stationary queue length distribution would in principle lead 
to a representation of the equilibrium distribution as an infinite linear combination of product 
forms. For an extensive treatment of the JSQ system using a generating function approach, the 
interested readers is referred to mm- An alternative approach that is not based on generating 
functions is the compensation approach [1]. This approach directly solves the equilibrium 
equations and leads to an explicit solution. The essence of the approach is to first characterize 
the product forms satisfying the equilibrium equations for states in the inner region of the 
two-dimensional state space and then to use the product forms in this set to construct a linear 
combination of product forms which also satisfies the boundary conditions. The construction is 
based on a compensation argument: after introducing the starting product form, new product 
forms are added so as to alternately compensate for the errors on the two boundaries. 

The compensation approach has been developed in a series of papers [BSIEIE] and aims at 
a direct solution for a class of two-dimensional random walks on the lattice of the first quadrant 
that obey the following conditions: 

(i) Step size: only transitions to neighboring states. 

(ii) Forbidden steps: no transitions from interior states to the North, North-East, and East. 

(hi) Homogeneity: the same transitions occur according to the same rates for all interior 
points, and similarly for all points on the horizontal boundary, and for all points on the 
vertical boundary. 

The approach exploits the fact that the balance equations in the interior of the quarter plane 
are satisfied by linear (finite or infinite) combinations of product forms, that need to be chosen 
such that the equilibrium equations on the boundaries are satisfied as well. As it turns out, 
this can be done by alternatingly compensating for the errors on the two boundaries, which 
eventually leads to an infinite series of product forms. The SED queueing system in this paper 
is more complicated than the JSQ and other classical queueing systems, see e.g. [BEllllISlE], 
since the two-dimensional random walk that describes the SED system exhibits inhomogeneous 
behavior in the interior of the quadrant. In this paper, we show that the compensation approach 
can nevertheless be further developed to overcome the obstacles caused by the inhomogeneous 
behavior of the random walk. This leads to a solution for the stationary distribution in the 
form of a tree of product forms. 

The only other work in this direction is [3], which considers SED routing for two identical 
parallel servers with Poisson arrivals and Erlang distributed service times. The crucial differ¬ 
ence with our setting is that we do not focus on generalizing service times, but instead consider 
servers with different service rates. 

The remainder of the paper is organized as follows. In Section [^ we introduce the model 
in detail and describe the equilibrium equations. We discuss the compensation approach and 
its methodological extensions together with our contribution in Section [^ Some numerical 
results are presented in Section [^ Section [^ applies the compensation approach to determine 
the equilibrium distribution of the SED system as a series of product-form solutions. Finally, 
we present some conclusions in Section]^ 
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Variable 

Expression 

Interpretation 

qi 


Number of customers (groups of size 1) in queue 1 

92 


Number of customers in queue 2 

[ 92 /sj 


Number of groups of size s in queue 2 

m 

m = min(gi, [ 92 /sj) 

Minimum number of groups in queue 1 and 2 

n 

n = [ 92 / 5 ] - 91 

Difference between number of groups in queue 2 and 1 

r 

r = mod(g' 2 ) s) 

Number of customers in queue 2 that are not in a group 


Table 1: Interpretation of the state variables of the Markov process. 


2 Equilibrium equations 


The Markov process associated with the SED system has an inhomogeneous behavior in the 
interior of the quadrant, specifically, along the line s{qi + 1) = gi + 1, see Figure In this 
section we transform the two-dimensional state space {qi, ( 72 ) to a half-plane with a finite third 
dimension. For this state description, we show that the theoretical framework of the compen¬ 
sation approach can be extended and in this way we determine the equilibrium distribution of 
the SED system. 

We will henceforth assume that s is a positive integer number. The service rate s could 
also be chosen to be rational and the analysis would be similar, but notationally more difficult. 


We further elaborate on this point in Remark 2.1 


In queue 2 we count the number of groups of size s and denote it as j, i.e. j = \_q 2 /s \, and 
we denote the number of remaining customers as r, i.e. r = mod(g 2 ; s). Clearly, a single group 
in queue 2 requires the same expected amount of work as a single customer in queue 1. The 
total number of customers in queue 2 is thus js + r and for an arriving customer the expected 
delay in queue 2 is j -|- (r -|- l)/s. In terms of these variables, SED routing works as follows: 
if (71 -|- 1 < j + {r + l)/s the arriving customer joins queue 1 and if -|- 1 > j -|- (r -|- l)/s 
the arriving customer joins queue 2. In case the expected delays in both queues are equal, 
i.e. -|- 1 = j -|- (r -|- l)/s, the arriving customer joins queue 1 with probability q and queue 2 
with probability 1 — q. 

For convenience we introduce the length of the shortest queue m = min(gi,j) and the 
difference between queue 2 and queue 1, i.e. n = j — qi. Using this notation, the SED system 
is formulated as a three-dimensional Markov process with state space {(m, n,r) | m G No, n G 
Z, r = 0,1,..., s — 1}. Under the stability condition (1.1) the equilibrium distribution exists. 


Let p{m, n, r) denote the equilibrium probability of being in state (m, n, r) and let 


p(m, n) = (p(m, n, 0),p{m, n, 1),... ,p(m, n,s — 1))^ 


( 2 . 1 ) 


with the transpose of a vector x. Throughout the paper, we use bold lowercase letters or 
numbers for vectors and uppercase Latin letters for matrices. For convenience, we have listed 
all state variables and their interpretation in Table 
The transition rates are given by 


(m, n, r) 


(l+s)p 


(m -I- l,n - l,r), 
< (m, n, r -|- 1), 

{m + l,n + 1,0), 


m > 0, 

n > 0, r 

m > 0, 

n < 0, r 

m > 0, 

n < 0, r 


(m, 0, s — 1) - — ^ (m, 1, 0), m > 0, 

(m, 0, s — 1) ——)■ (m, —1, s — 1), m > 0, 


0,1,... ,s - 1, 
0,l,...,s-2, 
s-1, 
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(a) Three-dimensional (m, n, r) transition rate di¬ 
agram. Note that the third dimension is perpen¬ 
dicular and gives rise to s-layers. 



(b) Partitioning of the state space of 
the Markov process: interior (or inner) 
states horizontal states =; and verti¬ 
cal states III. 


Figure 2: Transition rate diagram and state space partitioning. 


corresponding to arrivals, and 


m > 0, n > 0, r = 0,1,..., s — 1, 
m > 0, n < 0, r = 0,1,..., s — 1, 

m > 0, n £ Z, r = 1, 2,..., s — 1, 

m > 0, n > 0, r = 0, 

_ (m — 1, n — 1, s — 1), m > 0, n < 0, r = 0, 

corresponding to service completions. Figure [^a) displays the transition rate diagram for the 
three-dimensional state space. The transition rates are described by the matrices Ax^y in the 
positive quadrant and Bx^y in the negative quadrant, where the pair {x, y) indicates the step 
size in the (m, n)-direction. Let I be the s x s identity matrix, be an s x s binary 

matrix with element (x, y) equal to one and zeros elsewhere, and L an s x s subdiagonal matrix 

with elements {x,x — 1), x = 1, 2,..., s — 1 equal to one and zeros elsewhere. For consistency 
with the indexing of the vector p(m,n), indexing of a matrix starts at 0. The transition rate 
matrices take the form 


[m, n, r) 


[m, n, r) 




(m — 1, n -|- 1, r), 
(m, n -|- 1, r), 

(m, n,r — 1), 

(m, n — 1, s — 1), 


Ai-i = (1 + s)pl, 

^-1,1 = -Bo,i = I, 

^0,0 = ~(1 + 'S)(p + 1)-^ + sL^ 


Aop = (1 -h s)p(l - 
i?i,i = (l + s)pM(0’^-i), 

-^0,0 = ^0,0 + (1 + s)pL. 


The equilibrium equations can be written in matrix-vector form. We partition the state 
space as illustrated in Figure j^b). For the interior of the positive and negative quadrant we 
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have the following inner equations 


^o,oP(W', n) + 74i _ip(m - 1, n + 1) + Ao,_ip(m, n + 1) 

+74_i_ip(m + 1, n — 1) = 0 , m > 1, n>2, (2.2) 

n) + i?i 4 p(m - 1, n - 1) + i?o,iP(”i, n - 1) 

+i?_i^_ip(m + 1, n + 1) = 0, m > 1, n <—2. (2.3) 

The equilibrium equations corresponding to the states on the horizontal axis, or directly adja¬ 
cent to the horizontal axis, are referred to as the horizontal boundary equations and are given 
by 

1) -h _ip(m - 1 , 2 ) + Ao-ip(m, 2) 

+A_i^ip{m -M, 0) ^o,lP(”^, 0) = 0, 

^o,oP("2-, - 1 ) + Bi^ip{m - 1 , -2) So,ip(m, -2) 

+B-i-ip{m -I- 1, 0) -I- Bo^_ip{m, 0) = 0, 

Bofip{m, 0) _ip(m - 1,1) 5i,ip(m - 1, -1) 

+Ao-ip{m, 1) .Bo,ip(m, -1) = 0 , 

The vertical boundary equations are 

(^ 0,0 + -^)p(0, n) -I- ^o,-iP(0, n -I- 1) -I- ^_iqp(l, n - 1) = 0, m = 0, n > 2, (2.7) 

{Bofl + sM*^°’°^)p(0, n) -I- i?o,iP(0, n - 1) -|- i?_i,_ip(l, n -|- 1) = 0, m = 0, re < -2. (2.8) 

Finally, for the three remaining boundary states near the origin, we have 

(7lo,o + /)p(0, 1) + 7lo,-ip(0, 2) + ^-i,ip(l, 0) + 7lo,ip(0, 0) = 0, (2.9) 

(So,o + sM(0’0))p(0, -1) + 5o,ip(0, -2) + S_i,_ip(l, 0) + So,_ip(0,0) = 0, (2.10) 

{Bo,o + 1 + sM(°’0))p( 0, 0) + Ao,-ip(0,1) + Bo,ip(0, -1) = 0. (2.11) 


rre > 1, 

re = 1, 

(2.4) 

rre > 1, 

re = —1, 

(2.5) 

rre > 1, 

re = 0. 

(2.6) 


Remark 2.1. {Rational s) A system with a rational service rate s = ^ can also be analyzed. 
In that case, one needs to consider a system with two servers and service rates si and S2- 
Similar to our analysis at the start of Section one denotes the number of groups of size si in 
queue 1 as i and the number of groups of size S2 in queue 2 as j. Then, let G {0,1,..., — 1} 

denote the number of remaining customers in queue re = 1,2. Based on the aforementioned 
construction, a single group in either queue 1 or 2 requires the same expected amount of work. 
Lastly, set m = min(i, j), n = j — i and the third finite dimension is a lexicographical ordering 
of the states (ri, r 2 ) G {0,1,..., si — 1} x {0,1,..., S 2 — 1}. This state space description leads 
to a transition rate diagram that has a similar structure as the one seen in Figure [^a). In this 
sense, a system with a rational service rate can be analyzed using the approach described in 
this paper. 

3 Evolution of the compensation approach and our contribu¬ 
tion 

In this section we use the abbreviations: vertical boundary (VB); vertical compensation step 
(VCS); horizontal boundary (HB); and horizontal compensation step (HCS). 
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(a) Identical servers [4] (b) Non-identical servers (c) Identical servers and Er¬ 

lang arrivals with s phases 
(s-layered transition rate di¬ 
agram) [2] 

Figure 3: Simplified transition rate diagrams on the state space (m, n) for JSQ systems with 
two servers. 


The compensation approach is used for the direct determination of the equilibrium distri¬ 
bution of Markov processes that satisfy the three conditions mentioned in Section The key 
idea is a compensation procedure: the equilibrium distribution can be represented as a series 
of product-form solutions, which is generated term by term starting from an initial solution, 
such that each term compensates for the error introduced by its preceding term on one of the 
boundaries of the state space. In this section, we motivate why the SED system requires a 
fundamental extension of the compensation approach. We do so by first describing the evolu¬ 
tion of the compensation approach through a series of models and present for each model the 
corresponding methodological contribution. Finally, we describe the extension required for the 
SED system. 

The compensation approach was pioneered by Adan et al. [1], for a queueing system with 
two identical exponential servers, both with rate 1, and JSQ routing. Such a queueing system 
can be modeled as a Markov process with states {qi,q 2 ) £ Nq, where qi is the number of 
customers at queue i, including a customer possibly in service. By defining m = min(gi,g 2 ) 
and n = q 2 — qi, one transforms the state space from an inhomogeneous random walk in the 
quadrant to a random walk in the half plane that is homogeneous in each quadrant. Since the 
two quadrants are mirror images of each other, it is not needed to determine the equilibrium 
probabilities in both quadrants; it suffices to do so in the positive quadrant. The transition 
rate diagram of the Markov process is shown in Figure |^a). 

For the symmetric JSQ model in [3|, the initial solution is of the form riQa’^^Q, where % is 
a coefficient and ao and /3o are the compensation parameters that satisfy the kernel equation 

a/32{p+l) = ^'^2p + a^‘^+ a^, (3.1) 
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Figure 4: For a system with two identical servers and JSQ routing, the compensation approach 
generates in each compensation step a single compensation term. 


which is obtained by substituting a”*/?”' in the equilibrium equations of the interior of the 
positive quadrant and dividing by common powers. This initial solution satisfies the equilibrium 
equations in the interior and on the HB (there is only one such solution). In order to compensate 
for the error on the VB, one adds the compensation term such that r/o^cT/^o +^oct]"/^o 

satisfies the equilibrium equations in the interior and on the VB. The compensation parameter 
Oil with ai < /3o is generated from (3.1) for a fixed /3 = /3o- The coefficient uq satisfies a 
linear equation and is a function of r/o, ao, /5o) and ai. The resulting solution violates the 
equilibrium equations on the HB. Hence, one adds another compensation term such 

that r'o<a™/5o + satisfies the equilibrium equations in the interior and on the HB, 

where /3i and rji are determined in a similar way as for the VCS. Repeating the compensation 
steps leads to a series expression for the equilibrium probabilities that satisfies all equilibrium 
equations: 


OO 

p{m, n) = CY^ ^ X] rn > 0, n > 1, 


1=0 


(3.2) 


1=0 


where C is the normalization constant. Figure displays the way in which the compensation 
parameters are generated. 

The first extension is presented in [5] for the asymmetric JSQ model, i.e. the servers are 
now assumed to be non-identical with speeds 1 and s. The symmetry argument used earlier 
does not hold anymore and one needs to consider the complete half-plane, see Figure j^b) for 
the transition rate diagram. Note that the half-plane consists of two quadrants with different 
transition rates that are coupled on the horizontal axis. The approach in this case is an 
extension of the approach introduced in [5]. In a VCS, one compensates solutions that satisfy 
the positive inner equations on the positive VB as well as solutions that satisfy the negative 
inner equations on the negative VB. Two kernel equations (one for each quadrant) are used 
to generate the a’s, and the coefficients satisfy different linear equations. For a HCS, each 
product-form solution that satishes the positive inner equations, generates a single /3 for the 
positive quadrant and a single (5 for the negative quadrant. Accordingly, a product-form 
solution that satisfies the negative inner equations is compensated on the HB. Thus, in this 
case the generation of compensation parameters has a binary tree structure, see Figure 

A further extension of the compensation approach is presented in [2] for a model with 
two identical servers, Erlang-s arrivals and JSQ routing. The state description is enhanced by 
adding a finite third dimension that keeps track of the number of completed arrival phases. The 
random walk in the positive and negative quadrant are mirror images, which permits to perform 
the analysis only on the positive quadrant, see Figurej^c) for the transition rate diagram. In [2] 
the authors extend the compensation approach to a three-dimensional setting. Due to the three- 
dimensional state space, each compensation term takes the form a^f3'^i+{a, /3), where i+(a,/3) 
is a vector of coefficients of dimension s (equal to the number of arrival phases). In each HCS, 
s different parameters j3 are generated instead of just one. A graphical representation of the 
generation of compensation terms, which has an s-fold tree structure, is depicted in Figure 

Our contribution. The model at hand is defined on an s-layered half plane, thus requires 
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initial vertical horizontal 

solution compensation compensation 

Figure 5: For the asymmetric JSQ model, the compensation approach generates different 
compensation terms for the positive (straight arrow) and the negative quadrant (snaked arrow). 



initial vertical horizontal 

solution compensation compensation 

Figure 6; For the symmetric JSQ model with Erlang-s arrivals, the compensation approach 
generates s compensation terms in a HCS and just one compensation term in a VCS. 
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that we further extend the compensation approach. Similarly to [5], we need to account for 
the two quadrants by considering two kernel functions (one for each quadrant). Furthermore, 
in accordance with [2], in every HCS, a total of s different parameters j3 are generated for 
the positive quadrant and a single (3 for the negative quadrant. This leads to a (s + l)-fold 
tree structure for the compensation parameters as depicted in Figure Additionally, the 
product-form solutions take the form /3) or j3) depending on whether 

they are defined in the positive or negative quadrant, respectively. The resulting solution for 
the equilibrium distribution is, for m > 0, n > 1, 

oo (s-|-l)* s 

p{m,n) = CY^ 

1=0 i=l j=l 

oo (s+1)* s 

+ c EEE ^l+l,d{i)+j^l+l,d{i)+jf^l,d{i)+j^+(^l+'^3^ (^l,d{i)+j) ■ (3.3a) 

1=0 i=l j=l 

For m > 0, n = 0, 

oo (s-l-1)* 

p(m,n) = C'^ Y 

1=0 i=l 


For m > 0, n < —1, 


oo (s-l-1)* 

p(m,n) = C'^ Y 
1=0 i=l 
oo (s-l-1)* 

+ c EE i A,j(s-|-1)) 7 (3.3c) 

1=0 i=l 


where C is the normalization constant and d{i) := (i — l)(s + l). The first subscript I is the level 
at which a parameter resides, starting at level I = 0 (the initial solution). Within a level, the 
parameters are differentiated by using an additional index i. A horizontal compensation step 
and the initial solution has coefficients rj and a vertical compensation step has coefficients 
Additionally, a vector h is generated in each horizontal compensation step. The initial solution 
is described in Lemma 5.1H the horizontal compensation step is described in Lemma 5.13, and 
the vertical compensation step is described in Lemma 5.12, see Figure 


4 Numerical results 


Expression (3.3) is amenable for numerical calculations after applying truncation. For m > 

0, re > 1, 


1 2 j (s-l-1)* S 

PL{m,n) = CY 

1=0 i=l j=l 
L^J (s-Hi)' 

+^E EE ^l+l,d{i)+j^l+l,d{i)+jf^l,d(i)+j^+('^l-+'^d^ Pl,d{i)+j)- (4.1a) 

1=0 i=l j=l 
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initial vertical horizontal 

solution compensation compensation 


Figure 7: For the SED system, the compensation approach generates different compensation 
terms for the positive (straight arrow) or the negative quadrant (snaked arrow). The number 
of terms generated in the positive and negative quadrant are different. 
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Qi 


qi 


qi 


(a) p = 0.6 


(b) p = 0.75 


(c) p — 0.9 


Figure 8: Heat plot of the equilibrium distribution mass for the SED system with s = 3, 
q = 0.4, and varying p, determined according to (4.1) with L = 16. The heat plot shows 
where the probability mass is located (darker colour means more mass). The dashed lines are 
qi + l = {q 2 + l)/s and qi = q 2 /s. 


For m > 0, n = 0, 


If J (^+1)' 

PL{rn,n) = C'^ ^ (4.1b) 

1=0 i=l 


For m > 0, n < —1, 


LfJ (^+1)' 

PL(m,n) = C'^ ^ 

1=0 i=l 

+ ^ X] X] ^;+l,j(s+l)“I+l,i(s+l)/^Z,i”s+l)i-(“/+l,i(s+l)> A,i(s+1))) (4Tc) 

1=0 i=l 


where the empty sum Yl7=o Here, L = 0 indicates only the initial solution and for instance 
L = 3 indicates an initial solution, a vertical, horizontal and another vertical compensation. 
Naturally, as L increases, the approximation becomes more accurate. 

We perform several numerical experiments that verify that under SED routing the joint 
queue length process concentrates between the line where the expected delays in both queues 
are equal qi + 1 = {q 2 + l)/s and the line where the expected waiting time is equal qi = q 2 /s 
using the equilibrium distribution. To this end, we consider a system with service rates 1 and 
s = 3, q = 0.4 and set L = 16 in (4.1). The equilibrium distribution for this model and varying 
p is given in Figure]^ in the form of a heat plot, supporting our claim. 


Next, to demonstrate the rate of convergence of the series in ( |3.3[ ), we derive the number of 
compensation steps L for which the equilibrium probabilities p(m, n) are considered sufficiently 
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Figure 9; The number of compensation steps L for each state (m, n) such that the resulting 


PL{m,n) is accurate according to (4.2). The dashed line is m + \n\ = 3. Parameters are s = 4, 
p = 0.8, and q = 0.4. 


accurate: As a measure of accuracy we compute for each state (m, n) the minimum number of 
compensation steps L such that 


\pL{m,n,r) - pL-i{m,n,r)\ ^ .,„_4 

max -^-r- < 10 

r=o,i,...,s-i r) 


(4.2) 


Figure shows that away from the origin, the convergence of the series is very fast, but the 
convergence is also quite fast for states close to the origin. Note that the distance of a state 


(m, n) to the origin is directly related to the rate of convergence of the series expression (3.3). 
In particular, it seems to be a function of m + |n|: faster convergence further away from the 
origin. This property is formally proven in Section 5.7 and can be exploited for numerical 
computations. 

Remark 4.1. [No curse of dimensionality) Take a triangular set of states 7m = {{m,n) \ 
m G No, re G Z, rre + |re| < M}, where M is some non-negative integer. Technically, M needs 
to be strictly larger than some non-negative integer N, but we do not go into the details here; 


the lower bound N is described in Section 5.7 For states outside Tm, we can use (4.1) to 


compute p(rei, re). Since the number of compensation steps L required to achieve accurate 


results according to (4.2) decreases with m-|- |re|, the number of compensation steps L for each 
state [m,n) ^ Tm is relatively small. For states {m,n) G Tm, a linear system of equilibrium 
equations needs to be solved to determine the equilibrium probabilities p(?re, re), where one uses 
that p(rre, re), (m, re) ^ Tm are known. 

As an example, for s = 4, /? = 0.8 and q = 0.4 the choice M = 3 implies that L = 1, 
which gives accurate results for p(rre,re), [m,n) ^ 7m, see Figure]^ So it is evident that the 
compensation approach does not suffer from the curse of dimensionality. 
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5 Applying the compensation approach 


5.1 Outline 


In the following subsections we describe the main steps in constructing the equilibrium distri¬ 


bution of the SED system. First, we develop some preliminary results in Section 5.2, showing 
that the inner equations have a product-form solution and we determine one of the two param¬ 
eters of the product-form solution explicitly. Using these preliminary results, we determine the 


unique initial solution in Section 5.3 The vertical compensation step is outlined in Section 5.4 


Section 5.5 describes the horizontal compensation procedure. We formalize the resulting so¬ 
lution in terms of a sequence of product-forms in Section [5. 6[ This sequence of compensation 
terms grows as a (s -|- l)-fold tree and the problem is in showing that the sequence converges. 


Section 5.7 is devoted to the issue of convergence. 


5.2 Preliminary results 

We conjecture that the inner equations have a product-form solution. To this end, we examine 
a related model that has the same behavior in the interior as the original model. We then 
show that the equilibrium distribution of this related model can be expressed as a product- 
form solution. Moreover, modeling this process as a quasi-birth-and-death (QBD) queue, we 
obtain a closed form expression for one of the parameters of the product-form solution. This 
procedure is closely related to the one in [2l Section 4]. 

The related model is constructed as follows. We start from the state space of the original 
model in Figure [^a) and bend the vertical axis as shown in Figure 10 Note that for this 
modihed model, m G Z. We add an additional state (—1,0, s — 1) with a transition rate 
Ao,i = seo to state (—1,1) and a transition rate Bq-i = (1-|-to state (—1, —1), where 
e* is a column vector of zeros of length s with a one at position i. The transitions from the 
states on the diagonal m-|-n = 0, n>l are kept consistent with the transitions from a state in 
the positive interior of the SED model, where the downward transitions are redirected to the 
state (—1,0,5 — 1), specihcally, ^o,-i = sGq . Similarly, the transitions from the states on the 
diagonal m — n = 0, n < —1 are kept consistent with the ones in the negative interior of the 
SED model, where the upward transitions are redirected to the state (—1, 0, s — 1), specifically, 
Bq^i = 1^, where 1 is a column vector of ones of size s. 

The characteristic feature of the modified model is that its equilibrium equations for m + 
|n| = 0, |n| > 2 are exactly the same as the ones in the interior and the equilibrium equations 
for n G { — 1,0,1} are exactly the same as the ones on the horizontal boundary of the original 
model. In this sense, the modified model has no “vertical boundary” equations. 

We next present two lemmas. The first lemma states that a product-form solution exists 
for the modified model, while in the second lemma we identify the geometric term of the 
product form expression. Let p(m, n) = {p{m, n, 0),p{m, n, 1),... ,p{m, n,s — 1))^ denote the 
equilibrium distribution of the modified model. 

Lemma 5.1. For p < 1, the equilibrium distribution of the modified model exists and is of the 
form 

p(m,n) = a”^q(n), m-|-|n|>0, nGZ (5-1) 

with a G (0,1) and q(n) = {q{n, 0), q{n, 1),..., q{n, s — 1))^ such that 


OO 

E 


a 


r) < OO, r = 0,1,..., s — 1. 


(5.2) 
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Figure 10: s-layered transition rate diagram of the modified model. Note that the third 
dimension, i.e. r, is perpendicular to the page. The non-filled dot, graphed at position (—1,0), 
corresponds to the single state (—1, 0, s — 1). 


Proof. First notice that the modified model is stable whenever the original model is stable and 
that p(m, n) satishes the inner and horizontal boundary equations (2.2)-(2.6) for all m+|n| = 0, 
except for state (0,0). Observe that the modified model restricted to the area {{m,n) \ m G 
No, n G Z, m + |n| > no} U {(no — 1, 0, s — 1)}, no = 1, 2,... embarked by two lines parallel 
to the diagonal axes, yields the exact same process. Hence, we can conclude that 


p(m + 1, n) = Q!p(m, n), m -|- |n| > 0, n G Z 

and therefore 

p{m, n, r) = a"^q{n, r), m + |n| > 0, n G Z, r = 0,1,.. 
Finally, we observe that 


,s - 1. 


a l”lg(n,r)= p(—|n|, n, r) < 1, 


which concludes the proof. 


(5.3) 

(5.4) 

(5.5) 
□ 


In Lemma 5.1 we have shown that the equilibrium distribution of the modihed model has 
a product form which is unique up to a positive multiplicative constant. In the next lemma we 
determine the unique a in (5.1). More concretely, the unique a is equal to 

Lemma 5.2. For p < 1, the equilibrium distribution of the modified model is of the form 

p(m, n) =/9^^’'“^^™'q(n), m-|-|n|>0, n G Z. (5.6) 

Proof. Let k denote the total number of customers in the system, i.e. 


k = 


(1 -|- s)m + sn + r, n > 0, 
(1 -|- s)m — n + r, n < 0. 


(5.7) 
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Then, the number of customers in the system forms a Markov process and for all states k > s, 
the transitions are given by: (i) from state k to state k + 1 with rate (1 + s)p-, and (ii) from 
state fe + 1 to state k with rate 1 + s. Let pk denote the probability of having k customers in 
the system. From the balance principle between states k and /c + 1 we obtain 


Pk+i = ppk, k> s. 


(5.8) 


Furthermore, 


Pk-\-s-\-l — 


E 


p{m, n, r) + 


E 


{m,n,r) 

(H-s)mH-sn+r=fc+sH-l 


(m,n,r) 

{l-\-s)m—n+r=k+s-\-l 


E 


a^q{n, h) + 


E 


p{m, n, r) 
a'^q{n, r) 


{m,n,r) 

{l-\-r){m—l)-\-sn-\-r=k 


(m,n,r) 

(l-\-s)(m—l)—n-\-r=k 


= a a^q{n,h) + a a^q{n,r) 


{l,n,r) 

{l-\-s)l-\-sn^r=k 


{l-\-s)l—n-\-r=k 


= cxpk. 

Combining the last two results immediately yields a = 


(5.9) 

□ 


Combining Lemmas 5.1 and 5.2 yields the following result for the original model. 


Proposition 5.3. For p < 1, the inner and horizontal boundary equations (2.2)-(2.6) have a 
unique solution of the form 

p(m, n) =/3^^’''®^™'q(n), m > 0, nGZ (5.10) 

with q(n) = {q{n, 0), q{n, 1),..., q{n, s — l))'^ non-zero and 

OO 

p~^^~^^'^^"'^q{n,r) < OO, r = 0,1,..., s — 1. (5-11) 


The solution obtained in Proposition 5.3 satisfies the inner and horizontal boundary equa¬ 
tions. However, we still need to specify the form of the vector q(n). It will become apparent, 
from Lemmas 5.4 5.6 and 5.9 that the form of the vector q(n) is entirely different in the 
positive quadrant, on the horizontal axis, and the negative quadrant. Correctly identifying 
q(n) will result in the initial solution that satisfies the equilibrium equations of the interior 
and the horizontal axis. In the following lemmas we describe the form of a solution satisfying 
the inner equations. 

Lemma 5.4. 


(i) The product form p(m, n) = a™/3”i+(a,/3), 
equations of the positive quadrant (2.2) if 


m > 0 , n > 1 is a solution of the inner 


D+{a, /3)i+{a, /3) = 0 


(5.12) 


with 


P) — ®/3^o,o + cxP'^Aq—i A—ip -(- (5^A\—i 
and the eigenvector i+(Q!, (3) = (i+(a, /?, 0), z+(a, (3,1),, i+{a, /3,s — 1))"^ satisfies 

i+{a,f3,r) _ f af3{l + s){p + 1) - /3‘^{1 + s)p - 0 “^ y _nT 1 

i+{a,/3,0) V al3s j , r , ,...,s . 


(5.13) 


(5.14) 
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(ii) The product form p(m,n) = a™'/^ "'i_(a,/3), m > 0, n < —1 is a solution of the inner 
equations of the negative quadrant (|2.3|) if 


with 


D_{a, I3)i_{a, /3) = 0 

D_ {a, (3) = a/3-Bo ,0 + a/3^Bo,i + + 0^ Bi^i 


(5.15) 

(5.16) 


and the eigenvector i_(a, /3) = (i_(a, /3,0), i_{a, /3,1),..., i_{a, f3,s — 1))^ satisfies 
i_ (a, /3, r) ^(a, /3, V>_ (/3))V’+(/3)" - '&(«, /?, V'+(/3))V’- (/?)’' 


i_(a,/3,0) 

with 


V'±(/3) = 


^'(a,/3,V’-(/3)) - ^'(a,/3,V’+(/3)) 

(1 + s)(p + 1) - /3 ± ^(/3 - (1 + s)(/9 + 1))2 - 4s(l + s)p 


r = 0, l,...,s-1 (5.17) 


2s 


(5.18) 


and 


'^{a, f3,^p) = (3 - {I + s){p + 1) + sip +—{1 + s)p^0^ ^ = {l + s)pip ^(-'0^-1). (5.19) 

a a 


Proof, (i) Inserting the product form ip{m,n) = a™/3"'i_,.(a,/3) into 


and dividing by 


^m-ipn-i (5.12). For det(B+(a,/3)) = 0, the rank of D+{a,l3) is s — 1 and thus 

we have one free variable, which allows us to express i+{a,l3,r) in terms of i+(a,/3,0). The 
eigenvector \0a,l3) follows from the system of linear equations (5.12), which reads 

(a^ + /3^(l + s)/? —a/3(l + s)(p+l))z+(a, /3, r) + a/3si+(a, /3, r+ 1) = 0, r = 0,..., s — 2 (5.20) 


m—1 o—n-\-l 


and gives ( 5.14[ ). 

(ii) Inserting the product form p (m, n) = 0^(3 "'i_(a, j3) into ( |2.3[ ) and dividing by a™' ^j3 
results in (5.15). For det(B_(a,/3)) = 0, the rank of D_{a,l3) is s — 1 and thus we have one 
free variable, which allows us to express i_{a, /3, r) in terms of i_(a, /3, 0). The system of linear 
equations (5.15) reads 

(a/3^ — al3{l + s)(/9 + l))i_(a, /3,0) + a/3si_{a, [3, 1) + /3^(1 + s)pi_{a, f3, s — 1) = 0, (5.21) 

al3{l + s)pi_{a,(3,r - 1) + (a/3^ - a/3(l + s){p + l))i_{a, (3,r) 

+ a/3sz_(a,/3, r + 1) = 0, r = 1, 2,..., s — 2. (5.22) 

We construct a solution of the form i_{a,l3,r) = {c+ip+{l3Y + c_ip_{/3Y)i_{a, (3,0) and thus 
require c+ + c_ = 1. The two roots Y±{P) follow from substituting i_{a,l3,r) = V'’'i_(a,/3,0) 
in (5.22) and dividing by /3,0), yielding 

sY'^ + (/3 - (1 + s){p + 1))Y + (1 + s)p = 0. (5.23) 


The constants c+ and c_ follow from (5.21) and the simplification follows from the fact that 
i-{a,ft,r) = P,0) satisfies (5.23). □ 

Remark 5.5. {Normalized eigenvectors) Without loss of generality we henceforth set the first 
elements z+(a,/3,0) = i_(a,/3,0) = 1 for all a and pi. Since normalization of the equilibrium 
distribution follows afterwards, one can arbitrarily select the value of the first element of the 
eigenvectors i+(a,/3) and i_(a,/3). 
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We wish to determine the a’s and /3’s with 0 < |a|, \(3\ < 1, for which (5.12) and (5.15) 
have non-zero solutions i+(a,/3) and i_(a,/3), respectively. Equivalently, we determine a’s and 
/3’s for which det(Zl+(a,/3)) = 0 or det(Zl_(a,/3)) = 0. We first investigate the location and 
the number of zeros of det(Il+(a, /3)) = 0. 

Lemma 5.6. The equation det(Il+(a,/?)) = 0 assumes the form 


(a/3(l -|- s)(/9 -|- 1) — /3^(1 -|- s)p — a^Y — P{aj3sy = 0 
and can he rewritten to 

a/3il +s){p + l) - Y^{1 +s)p-a'^ _ 
aps 

where Ui is the i-th root of unity of = 1 and is the principal root. 


(5.24) 


(5.25) 


(i) For every a with |a| G (0,1), equation (5.25) has exactly one root fdi inside the open circle 
of radius |a| for each i. Furthermore, all s roots /3j are distinct. 


(ii) For every ft with |/3| G (0,1), equation (5.25) has exactly one root a* inside the open circle 
of radius |/3| for each i. Furthermore, all s roots a* are distinct. 


Proof. Dividing (5.24) by (aftsY and taking the s-th root reduces (5.24) to (5.25). 
(i) The proof consists of three main steps: 


(a) For every fixed a with |a| G (0,1), equation (5.24) has exactly s roots (3 inside the open 
circle of radius |a|. 

(b) These s roots j3 are distinct. 


(c) For every fixed a with |a| G (0,1), equation (5.25) has at least one root /3j inside the 
open circle of radius |a| for each i. 


Combining these three steps proves that for every fixed a with |a| G (0,1), equation (5.25) has 
at exactly one root Yi inside the open circle of radius |a| for each i and the /3j’s are distinct. 
(i)(a) Equation (5.24) is a polynomial of degree 2s in /? and we will show that exactly s roots 
are inside the open circle of radius |a|. Other possible roots inside the open unit circle appear 
not to be useful, since these will produce a divergent solution, as follows from (5.11). Divide 
both sides of (5.24) by and set 2 ; = /3/a to obtain 


f{zY - as^z^+^ = 0 

with f(z) = (1 -|- s)(/9 -|- l)z — (1 + s)pz‘^ — 1. Then f(z) has the two roots 

p + i± P{p+ly - 4p/(i s) 


v± = 


2p 


(5.26) 


(5.27) 


Observe that, for 0 < p < 1, > 1 and 0 < < 1. Furthermore, for \z\ = 1 one establishes 

\fiz)Y = 1(1 + s)ip + l)z - (1 + s)pz^ - ir 

> 1(1-1-s)(p-|-l)| 2 ;| — (1-1- s)p\z\^ — 1|'^ = s"^ > |a|s®. (5.28) 

Hence, by Rouche’s theorem (see e.g. [25l Theorem 9.3.2]), equation (5.26) has exactly s roots 


inside the open unit circle. Thus, (5.24) has for each fixed |a| G (0,1) exactly s roots /3 inside 
the open circle of radius |a|. 
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(i)(b) Label the left-hand side of (5.26) as g{z). For a root to be of (at least) multiplicity two, 
it must hold that g{z) = g'{z) = 0, which gives 


{p + l)z -h (s - l)pz^ -1 = 0 


(5.29) 


with solutions 


= 


-{p -H) ± V(/9 -F If 4(s - l)p 
2{s-l)p 


(5.30) 


Note that for 0 < /o < 1, a;_ < 0, 0 < a;+ < 1 and 
of p with a;+ = 1 for /9 = 0 and a;+ = {yfs — l)/(s — 1) 
reveals the contradiction 


is a monotonically decreasing function 
for p = 1. For z = uj+ equation (5.26) 


_ {l + p{l-2u:^)Y 

0J+ 


(5.31) 


Thus, the s roots j3 inside the open circle of radius |a| are distinct. 

(i) (c) This step is presented in Appendix [A| The proof was communicated to us by A.J.E.M. 
Janssen. 

(ii) Set z = ajYi multiply (5.25[) by zs and rearrange to obtain 


z{{l + s){p + 1) - Ui/3^^''s) = z^ + {l + s)p. 


(5.32) 


Label the left-hand side as f{z) and the right-hand side as g{z). Note that f{z) has a single 
root within the unit circle. For |z| = 1 one establishes 

\f{z)\ = | 2:((1 -h s)(/5 + 1) - = 1(1 -h s)(p-M) - 

> (1 s)(p-F 1) - = (1 s)p-M -b s(l - 

> (1 + s)p -|- 1 > \z'^ -|- (1 -|- s)/9| = |5'(^)|. (5.33) 


Hence, by Rouche’s theorem we establish that (5.25) has exactly one root a* inside the open 
circle of radius 


for each i. These roots a* are distinct, since the Ui are distinct. 


□ 


Remark 5.7. {Identical eigenvectors) We note that for fixed 13 with |/3| G (0,1) and fixed i, 
(5.25) is a quadratic equation in a and has two solutions, say a+ and a_, satisfying |a_| < 
|/3| < |q;+|. Then, combining (5.14) and (5.25) we have the property that i+(a+,/3) = i+(a_,/3), 
which we will use later. 


We next provide information on the location and the number of zeros of det(Zl_(Q!, f3)) = 0- 
The polynomial form of det(Il_(a,/?)) = 0 is too complicated for application of Rouche’s 
theorem. Therefore we apply the following matrix variant of Rouche’s theorem. 

Theorem 5.8. (de Smit [H]) Let A{z) = {aip{z)) and B{z) = {bij{z)) be complex n x n 
matrices, where B{z) is diagonal. The elements aij{z) and bij{z), 0<i,j<n— 1 are 
meromorphic functions in a simply connected region S in which T is the set of all poles of 
these functions. C is a rectifiable closed Jordan curve in S \ T. Let xb and xa+b be the 
number of zeros inside C of det{B{z)) and det{A{z) + B {z)), respectively, andys andpA+B the 
number of poles inside C of det{B{z)) and det(A( 2 ;) -|- B{z)) (zeros and poles of higher order 
are counted according to this order). If \bi^i{z)\ > ^ alH = 0, 1,..., n — 1, 

then 

xa+b - VA+B =XB-yB- (5.34) 
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Lemma 5.9. 

(i) For every a with |a| G (0,1), the equation det(Zl_(a,/3)) = 0 can he rewritten to 

+ /3^((1 + s)py - (/?)"* + V’-(/3)*) = 0, (5.35) 


where iIj±{( 3) is given in (5.18), and has exactly one root fi in the open circle of radius 

|a|. 


(ii) For every (3 with |/3| G (0,1), equation (5.35) has exactly one root a in the open circle of 
radius |/3|. 

Proof, (i) Setting z = /3/a, we observe that 


D_{a, /3) = a‘^{A{z) + B{z)) 


(5.36) 


with 


A(z) = (1 + s)pzL + szL'^ + (1 + s)pz^M^^'^-^'^ + (5.37) 

B{z) = z{az — (1 + s)(/9 + 1))/. (5.38) 

Hence, det(Il_(a, /3)) = det(H( 2 ;) + H(z)). Let us take C to be the unit circle. Furthermore, 

we can verify that, for |a| G (0,1) and \z\ = 1, 

s-l 

= 1(1 + s){p + 1) - azi > (1 + s)/9 + s = ^ \aij{z)\, z = 0,1,..., s - 1 (5.39) 

j=0 

and the number of zeros and poles inside C of det(il( 2 ;)) is xb = s and ys = 0, respectively. 


The number of poles of det(H(z) + B{z)) inside C is pa+b = 0. By Theorem 5.8 we derive that 
det(H(z) + B{z)) has exactly s roots inside C. Moreover, 

det(H( 2 :) + B{z)) = {-zy~^ + z‘^{{l + s)p)® 

L./2J 


— z 


E(-F 


i=0 


S — t 


S — I 


{s{l + s)py{{l + s){p + l)- azy-^^ ). (5.40) 


So, det(H(z) + B{z)) has the root 0 of multiplicity s — 1 and exactly one non-zero root inside 


C. The summation in (5.40) is known as the Waring formula. By using the following identity: 

L./2J 




i=0 


s — I 


s — l 


{xyy{x + y) 


s—2i 


(5.41) 


we can simplify the Waring formula, see e.g. 1201 equation (1)]. So, we set xy = s(l -|- s)p and 


X + y = [1 + s){p + 1) — cxz. Recall from (5.23) that ip+{az)'ijj_{az) = (1 -|- s)p/s. This allows 
us to establish that x = sy+{az) and y = sy_{az). Thus, (5.40[) simplifies to 


det(H(z) -|- B{z)) = {—zy ^ -|- z^((l -|- s)py — zs^(?/)+(az)^ -|- y_{azy'^ . (5.42) 

Since det(H(z) -|- B{z)) has exactly s roots in the unit circle and 0 is a root of multiplicity 
s — 1 we are led to the conclusion that there exists exactly one non-zero root in the unit circle. 


Thus, (5.35) has exactly one non-zero root (3 in the open circle of radius |a|. 


(ii) The proof is identical to the proof of (i). 


□ 


Remark 5.10. {Notation) For a fixed a with |a| G (0,1) the roots of (|5.24 ) are denoted by 
(3i,I32, ■ ■ ■, Ps with |/3i| < |a|, i = 1,2,... ,s and the single root of (5.35) is denoted by I3s+i 


with |/3,j+i| < |a|. Similar indexing exists for the roots a for a fixed (3 with \j3\ G (0,1). 
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5.3 Initial solution 


Lemmas 5.4, 5.6 and 5.9 characterize basic solutions satisfying the inner equations (2.2)-(2.3). 
In Lemma 5.11| we specify the form of q(n). 

Lemma 5.11. (Initial solution) For |a| G (0,1), let /3i,/32, ■ ■ ■,/?« be the roots of (5.24) with 
IPil < |a|, i = 1,2, ...,s and i+{a,f3i) the corresponding non-zero vectors satisfying (5.14). 
Symmetrically, let (dg+i be the root of (5.35) with |/3s+i| < |a| andi_{a, fds+i) the corresponding 
non-zero vector satisfying (5.17). Then there exists exactly one a with |a| G (0,1) for which 
there exists a non-zero vector h = (/i(0), h{l ),..., h{s — 1))^ and coefficients r]i,r] 2 , ■ ■ ■, Vs+i 
such that 




m > 0, 

n > 1, 


(5. 

II 

a^h, 

m > 0, 

n = 0, 

(5.43) 



m > 0, 

n < —1, 



satisfies (2.2)-(2.6). This unique a is equal to . 

The vector h is given by 

S 

h = a(^o,i + (dj) 

i=l 

and the coefficients r]i,r} 2 , ■ ■ ■ ,ils satisfy 

S 

+ uAq^i) + q;^So,o(^o,i + ^)i+(a,/3j) 


i=l 


= + aBo,i)i_(a,/3s+i) 


(5.44) 


(5.45) 


and T/s+i = 1. 


Proof. Substituting (5.43) into (2.4) yields 


1 * 

h =-(^0,1 + aA-i^y^ -i + af3fAo-i)i+{a,/di). (5.46) 

“ i=i 


Furthermore, by using (5.15) for the summands one obtains (5.44). Equation (5.45) follows 


from substituting (5.43) into (2.6) and using (5.44). One can arbitrarily choose rjs+i due to 
the normalization that follows at the end. □ 


5.4 Compensation on the vertical boundary 

The initial solution ( 5.43| ) does not satisfy the positive and negative vertical boundary. We 
consider a term r/a™’/3'^i+(a,/3) that satisfies the positive inner balance equations and stems 
from a solution that satishes both the inner and horizontal boundary equations. For now, we 


refer to this term as the original term. In case of the initial solution (5.43), this would be one 


of the s product-form solutions of the positive quadrant. The idea behind the compensation 
approach is to add compensation terms ^de original term, such that 

the linear combination 


rja 


'■/3''i+(a,/3)-b m > 0, n > 1, 


(5.47) 


2=1 
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satisfies both (2.7) and (2.2). Substituting (5.47) in (2.7) gives 


S 

r?/3'"“V+(a,/3)i+(a,^) + V V+(ai,/3)i+(ai, = 0 , n > 2, 


i=l 


(5.48) 


where V+(a,/3) = /3(Ao,o + 7) + /3^^o,-i + Note that we indeed require that the 

original term and the compensation terms share the same (5 and therefore we obtain the s 
roots ai, Q! 2 , ..., with |Q!i| < |/3| from (5.24). Clearly, (5.47) now satisfies the positive inner 
equations (2.2). This leaves the s coefficients Vi to satisfy the s equations (5.48). However, 


from Remark 5.7 we deduce that there exists a specific i for which i+(a,/3) = i+(ai,/3) and 
thus there is only one coefficient that is required to be non-zero. 

Let us now consider a term from the negative quadrant /3) that satisfies the 

negative inner balance equations and stems from the same solution that satisfies both the inner 
and horizontal boundary equations. Let us now refer to this term as the original term. In case 


of the initial solution (5.43), this would be the solution on the negative quadrant. As before, 
the compensation approach dictates to add a compensation term fi) to 

the original term, such that the linear combination 


r]a^l3 ”i_(a,;0)-h z^s+ia^i/3 "i_(Q;s+i, ;0), m > 0, n <-1, 
satisfies both (2.8) and (2.3). Substituting ( |5.49 ) in (2.8) gives 

(3)\_{a, j3) + («*+!, ,0) = 0, n < -2, 


(5.49) 


(5.50) 


where R-(a, (3) = -|- /3^Ro,i + Note that we indeed require that the 

original term and the compensation term share the same /3 and therefore we obtain the root 
Os+i with Icts+il < \/3\ from (5.35). This ensures that the linear combination (5.49) satisfies the 


negative inner equations (2.3). Even though there is only one coefficient Vs+i for s equations, 


it turns out that this provides enough freedom to satisfy (5.50) 
Lemma 5.12. (Vertical compensation) 


(i) Consider the product form f3) that satisfies the positive inner equations (2.2) 


that stems from a solution that satisfies the inner and horizontal boundary equations. For 
this a and fixed fi, let ai he the root that satisfies i+(a,/3) = i+(ai,/3) with |ai| < |/3|. 
Then there exists a coefficient ni such that 


p{m,n) = r]a'^f3'^i+{a, fi) + nia'fi'l3"'i+{ai, f3), m > 0, n > 1, 


satisfies (2.2) and (2.7). The coefficient ni satisfies 

i-g(i + »V 


= - 7 ] 


(5.51) 


(5.52) 


11 


Consider a product form (5) that satisfies the negative inner equations (2.3) 

that stems from a solution that satisfies the inner and horizontal boundary equations. 
For this fi, let Us+i he the root of (5.35) with |as+i| < \(3\ and let i_(as+i,/3) he the 


corresponding non-zero vector satisfying (5.17). Then there exists a coefficient i/g+i such 
that 


p{m,n) = ”'i_(a,/3) + fi), m > 0, n <-1, (5.53) 


22 












































satisfies (2.3) and (2.8). The coefficient is given by 


^ “ 1(1 + s)pi_{a,ffis - 1) 

Us + l = —11 -o-• 

® “ dTr(^ + s)pi-{as+iffi, s-1) 


(5.54) 


Proof, (i) Use (5.12) to establish that 


V+{a, I3)i+{a, fi) = {fil - —Ai_i)i+{a,/3) = fi{l- -{I + s)p)i+{a, fi). (5.55) 

a a 


Then using Remark 5.7 we can establish that there exists a single Vi that is non-zero. We label 


the single non-zero coefficient as vi and hnd it from (5.48) with the help of (5.55). 


(ii) Use (5.15) to establish that 

U_(a, fi) = — — (1 -|- s)/9M^°’^“^^)i_(a, j3). 


Thus, (5.50) reduces to a single equation and Ug+i follows directly. 


(5.56) 

□ 


5.5 Compensation on the horizontal boundary 


The solution obtained after compensation on the vertical boundary, as outlined in Lemma 5.12 


does not satisfy the horizontal boundary equations. So, we need to compensate for the error 
on the horizontal boundary by adding new terms. The compensation procedure on the hor¬ 
izontal boundary has a few differences from the one described in the previous section. We 
outline these differences by informally treating the compensation of a product-form term of 
the positive quadrant. The difference in the compensation procedure is due to the fact that 
adding compensation terms only for the positive quadrant does not make the solution satisfy 
the horizontal boundary equations. Intuitively, this originates from the fact that the horizontal 
boundary, i.e. n = 0, is connected to both the positive and negative quadrant. Thus, for the 
horizontal compensation step, we need to add product-form terms for the complete positive 
half-plane. It turns out, that these product-form terms are nearly identical to the initial so¬ 
lution, which indeed satished both the inner and horizontal boundary equations. The same 
procedure holds for a product-form term of the negative quadrant. 

The formal compensation on the horizontal boundary is outlined in the next lemma. 

Lemma 5.13. (Horizontal compensation) 


(i) 


Consider the product form that satisfies the positive inner equations (2.2) 

that stems from a solution that satisfies the inner and vertical boundary equations. For 
this a, let ffi, ■ ■ ■, ffi be the roots of (5.24) with \fii\ < |q;|, i = l,2,...,s and let 
i+{a,l3i) be the corresponding non-zero vectors satisfying (5.14). Symmetrically, for this 
a, let /3s+i be the root of (5.35) with |/3s+i| < |a| and let i_{a, (3s-{-i) be the corresponding 
non-zero vector satisfying (5.17). Then there exists a non-zero vector h and coefficients 
r]i,r] 2 , ■ ■ ■, Ps-i-i such that 


i/a™'/?”i+(a,/3)r/ia”^/3fi+(a,/3i), m > 0, n > 1, 
p(m,n) = a^h., m > 0, n = 0, 

m>0, n<-l. 


(5.57) 


satisfies (2.2)-(2.6) 
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The vector h and the coefficients ryi, r/ 2 , • • •, rjs+i satisfy 


(^0,1 + a^-i,i)h - a'^r]ii+{a, Pi) = z/m+(a,/ 3 ), 


(5.58) 


i=l 


aBoflh + rjiPi (^i-i + a^o-i)i+(a, A) 


1=1 


+r]s+iPs+i{Bi^i + aBo^i)i_{a, Ps+i) = -i + a^o _i)i+(a,/5), (5.59) 

— r/s+ias + ashifS) + (1 + s)pqh{s — 1) = 0. (5.60) 

(ii) Consider the product form P) that satisfies the negative inner equations 


(2.2) that stems from a solution that satisfies the inner and vertical boundary equations. 


For this a, let Pi,P2, ■ ■ ■, Ps be the roots of (5.24) with |/3j| < |a|, z = 1, 2,..., s and let 


i+{a,Pi) be the corresponding non-zero vectors satisfying (5.14). Symmetrically, for this 


a, let Ps-\-i be the root of (5.35) with |/?<j+i| < |a| and let i^{a, Ps+i) be the corresponding 


non-zero vector satisfying (5.17). Then there exists a non-zero vector h and coefficients 
r]i,r] 2 , ■ ■ ■, such that 


m > 0, n > 1, 

p(m,n) = m > 0, n = 0, 

va'^P~'^i_{a,P)+ r]s+ioT"Pf+P-.{a,Ps+i)-, m > 0, n <-1, 


(5.61) 


satisfies (2.2)-(2.6). 


The vector h and the coefficients r]i,r} 2 , ■ ■ ■, Vs-i-i satisfy 

S 

(^0,1 + a^-i,i)h - a'^r]ii+{a, Pi) = 0, 


(5.62) 


2 = 1 


aBQflh + ^ piPi {Ai_i _i) i+ {a, Pi) 


2=1 


+ris+iPs+i[Bi^i +aBo^i)i_{a,Ps+i) = -izP{Bi^i + aBo^i)i_{a, P), 
—ps^ias + ash{0) -|- (1 -|- s)pqh{s — 1) = uas. 


(5.63) 

(5.64) 


Proof, (i) Equations (5.58) and (5.59) follow from the exact same reasoning as outlined in 


the proof of Lemma 5.11 Equation (5.60) follows from substituting (5.57) in the negative 


horizontal boundary equations (2.5) and using (5.15). Note that this indeed reduces to a single 
equation. 

(ii) The proof is identical to the proof of (i). □ 


5.6 Constructing the equilibrium distribution 

The compensation approach ultimately leads to the following expressions for the equilibrium 
distribution of the SED system, where the symbol oc indicates proportionality. 
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(i) For m > 0, n > 1, 


OO {'5 + 1)* S 

p(m, n) oc E E 

1=0 1=1 j=l 

OO (s+1)^ s 

+ EEE i^i+i,d(i)+i«;+i,(i(i)+jA"d(i)+ji+(“«+i,i) Pl,d(i)+j)- (5.65a) 

/=0 i=l j=l 


(ii) For m > 0, 


OO (5 + 1)* 

p(m,0) oc ^ ^ 
1=0 i=l 


h 


l,i ■ 


(5.65b) 


(iii) For m > 0, n < — 1, 


OO (5+1)* 

p(m, n) oc EE i?z,i(s+i)«™i/5;_j(^+i)i-(a«,i, A,i(s+i)) 

Z=0 i=l 

OO (s+1)* 

+ EE l^Z+l,i(s+l)Clz+l,i(s+l)/^;^j(s_|_l)i- (ciZ+l,i(s+l); Pl,i{s+1))- (5.65c) 

/=0 i=l 


We briefly describe the indexing of the compensation terms, which grows as a tree. We indicate 
the level at which a parameter resides with I, starting at level I = 0. Within a level, we 
differentiate between parameters by using an additional index i. The procedure for generating 
terms is as follows: 


(I) The initial solution is determined from Lemma 5.11 


(V) For a vertical compensation step with fixed 

(i) If the index i is not a multiple of s +1, then the compensation terms are determined 
according to Lemma 5.12| (i) . 

(ii) Otherwise, the compensation terms are determined according to Lemma 5.12K ii). 
(H) For a horizontal compensation step with fixed 

(i) If the index i is not a multiple of s +1, then the compensation terms are determined 


according to Lemma 5.13 


(ii) Otherwise, the compensation terms are determined according to Lemma 5.13| (ii). 

Recall the definition d{i) := {i — l)(s + 1) and note that d{i) + s + 1 = f(s + 1). Informally, 
the indexing of subsequent terms is visualized in Figure [H Note that the first term in the 
compensation approach is ao,i = c.f. Lemma 


5.11 
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horizontal 

compensation 


vertical 

compensation 


horizontal 

compensation 


Figure 11: Indexing of the compensation terms for a subtree. A dashed rectangle indicates 
according to which lemma the compensation terms are generated. A straight arrow indicates 


that the roots a or /? are obtained through (5.24) and a snaked arrow through (5.35). 
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5.7 Absolute convergence 

In this subsection we prove that the series for the equilibrium probabilities, as formulated 


in (5.65), are absolutely convergent. Before we are able to do that, we need some pre¬ 


liminary results. Specifically, we establish that the sequences {az,i}igNo, i=i, 2 ,...,(s+i)' 
{A,i}zeNo, i=i, 2 ,...,(s+i)*+i decrease exponentially fast and uniformly in the levels of the param¬ 
eter tree. Furthermore, we investigate the asymptotic behavior of (ratios of) the parameters 
the coefficients, and the elements of the (eigen)vectors. 

Corollary 5.14. Let ai := A •= 2 ,...,(s+i)'+i IA,j| • Then, 

(i) = |ao,i| > /3o> ai> Pi > 012 > f32> ■■■ 

(ii) There exists c £ (0,1), such that 0 < ai,Pi < P . 


Proof, (i) Follows immediately from Lemmas 5.6 and 


(ii) For a fixed a, let Pi, P 2 ,..., /3s be the roots of (5.24) with |/3j| < |a|, i = 1,2,..., s and let 


5.9 


Ps+i be the root of (5.35) with |/3s+i| < |a|. We define 


t(a) := max \Pi/ot\ 

4 = 1,2,...,s + l 

— /ii+s 


(5.66) 


and let C = {a G C I |a| < ao,i = Using Lemmas 5.6 and 5.9 we have t{a) < 1 and 


in particular, since p £ ( 0 , 1 ), we have that C is a closed proper subset of the unit disc and 
thus Cl := maxcgc t(cK) < 1 - On® can reason along the same lines and obtain a bound C 2 for 
a fixed /3. The exponential decrease in terms of the level I follows from the fact that each /3 
that is generated from an a satisfies |/3| < |a|ci and each a that is generated from a /3 satishes 
|a| < \P\c 2 . Thus, we define c := C 1 C 2 . □ 


In Corollary 


5.14 


we have established that both sequences i=i, 2 ,...,(s+i)' ^cid 

{Pi,i}i£No, 4 = 1 , 2 ,...,(s+i)*+i tend to zero as / —>• 00 . Thus, letting a 4 , 0 or /3 0 is equiva¬ 

lent to letting / — 7 - 00 . In what follows, whenever the specific dependence on the level I is not 
needed, we opt to use the simpler notation of Sections |5.2H5.5[ The next lemma presents the 
limiting behavior of the compensation parameters and their associated eigenvectors. 

Lemma 5.15. 


(i) For a fixed a, let Pi,P 2 , ■ ■ ■, Ps be the roots of (5.24) with \Pi\ < |a|, i = 1 , 2 ,... ,s and 


let Ps+i be the root of (5.35) with |/3s+i| < |a|. Then, as a fO, 


(a) The ratio Pi/a -£■ v_, i = 1,2,..., s with v_ < 1, where v_ is the root of 

i;2(l-hs)p-i;(l-bs)(p-M)-hl = 0. (5.67) 


(b) The eigenvector i+(a, Pi) -£■ gq , i = 1,2,..., s . 

(c) The ratio Ps+i/a — )■ with < 1, where is the root of 

w;2((l + s)p)* - ws^{^P+{0)^ + V’-(O)*) + s'* = 0, 


(5.68) 


where V’±(0) is defined in (5.18). 


(d) The elements of the eigenvector i_(a, Ps+i,r) —i//g(0)^, r = 0,1,..., s — 1. 


(ii) For a fixed P, let ai,a 2 , ■ ■ ■ ,as be the roots of ( |5.24 ) with \ai\ < |;0|, i = 1,2,..., s and 
let cis+i be the root of (5.35) with |q;s+i| < |/3|. Then, as P ft), 
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(a) The ratio ai/P —)> l/f+, i = 1,2 ..., s with v+ > 1, where v+ is the root of (5.67). 

(b) The eigenvector /3) —>- 60 , i = 1,2 ..., s. 

(c) The ratio as+i/P —)■ l/w+ with w+ > 1, where w+ is the root of (5.68). 

(d) The elements of the eigenvector i_{as+i, (3, r) —>■ 'tp_{Oy, r = 0,1,..., s — 1. 


Proof. (i)(a) Set v = (3i/a and let a 0 in (5.26) to establish ( 5.67| ). The roots of (5.67) are 
dehned in (5.27) and satisfy 0 < < 1 and v+ > 1. 

(i)(b) We have that i+{a,/3,0) = 1 for all a and for a 4 , 0 the other elements in (5.14) go to 
((1 + s )(/9 + 1 ) — v_{l + s)p — l/v_)/s, which is equal to zero by (i)(a). 

(i)(c) Set w = j5s+i/a, divide ( 5.35| ) by and let a | 0 to establish (5.68). We note that if±{-) 
is a continuous function, so that indeed for a 4 , 0, if±{aw) V’±(0)- Furthermore, (5.68) is a 
quadratic equation in w with roots |rc_| < 1 and |u)+| > 1, where the inequalities follow from 
Lemma 15.91 and also satisfy w_,w+ G M+. 


(i)(d) Again, set w = fls+i/a. Recall the dehnition of the eigenvector in (5.17). We aim 
at showing that for a | 0, 4 '(q;,/ 3s+i, V^-(arc)) —)■ 0 and 4 '(q;,/ 3s+i, V^+(arc)) tends to some 
non-zero constant. As a 4 . 0, we have that 


^{a, /3s+i, Ip) ^ {1 + s)py ^{w.iy-l). 


(5.69) 


Note that V'_(0) G (0,1) and fj+iO) > 1. Since w_ < 1, the limiting value in (5.69) can only 
equal zero in the case that = l/'ip+{0Y. Thus, we have established that ^{a, fts+ijip-iaw)) 
tends to a non-zero constant as a 4- 0. We next verify that w_ = l/V’+(0)® is a solution to 
(5.68), which proves that T(a,/3s+i,'0+(at(;)) —)■ 0 as a 4 , 0. Substituting w = l/^p+{0)‘^ in the 
left-hand side of (5.68) yields 

:s%ip4or + y.{oY) + sy 


1 

V’+(o) 


2s ~ 


1 


1 


1 


^^+( 0 ) 


2s 


(((1 + s)pr - s^y^oYiMoy + Y-m + 

[iii+.s)pr-sYY+mY-io)r)=o, 


where we used V'+(O)'0-(O) = (1 -|- s)p/s, as found in (5.23). 
(ii) The proof is identical to the proof of (i). 


(5.70) 


□ 


Finally, we describe the limiting behavior of the coefficients. In the following lemma we 
introduce the variable 7 associated with a product-form solution that satishes the positive 
inner equations, and the variable 9 associated with a product-form solution that satishes the 
negative inner equations. 

Lemma 5.16. 

(i) Consider the setting of Lemma 5. M i). Then, os /? 4- 0, 

ni l-v_{l + s)p 


=■ lu¬ 


ll 


T] 1 - U+(l -I- s )/3 
Consider the setting of Lemma 5. T^ii) . Then, os /? 4- 0, 
Vs+i s -w_{l + s)pY+{0Y~^ 


(5.71) 


V 


s - t(;+(l -7 s)pY-{0) 


s-l 




(5.72) 
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Consider the setting of Lemma 5.13[i) 

Then, os a | 0, 

(a) 

h 

—y 0 


u 

(b) 

Vs+l ^ 

v_ — v+ 

v_s 

^ + w_y+{oy-^ 

(c) 



lim 


7 ??s + l- 


lojWI =: 7r?, 


(5.73) 

(5.74) 

(5.75) 


where Aj = (aj( 0 ), aj(l),..., aj(s —l))'^ is the solution to the following linear system 
of equations for a speeific j, with bj a eolumn vector of size s with unknowns, 


.Wa.j + Lhj = -v+Wj - ^eo 


-Wsij + (1 + s)p(l - = 


W4 


(5.76) 

(5.77) 


with 


( 1 


W := 


1 \ 


l/s 

V_ Ul 


1/s 

V_ U 2 


\v_ 


(s-l)/s s_l (s-l)/s s_l 


Ul V_ 


Ur, 


1/s 

V_ Us 


{s — l)/s Q _1 

vl " u! ^ 


w,' := ( 1 v^l^Uj 


{s—l)/s s-1 

V\ Uj 


(5.78) 


(5.79) 


(iv) Consider the setting of Lemma 5.13K ii). Then, as a | 0, 
(a) 


h 

-^ 0 . 

V 


(b) 


T]s+1 V.s^ +W+'lp_{ 0 y 


f_sbj + W_'lp+{ 0 )' 


=: a 


• ''Vs + l- 


(c) 


lim 

04.0 


< max |c(r)| =; 

0<r<s-l ' 


(5.80) 

(5.81) 

(5.82) 


where c = (c( 0 ), c(l),..., c(s — 1 ))^ is the solution to, with d a column vector of size 
s with unknowns. 


v_Wc + Ld = —{/w. 
-Wc + (1 + s)p(l - g)M(°’"-bd = 0. 


.y-(oy ^ + dn^+iw-y+io) 


^)eo, (5.83) 
(5.84) 


Proof, (i) Divide both sides of (5.52) by rj and let a | 0. 


1 . 


(ii) The proof is identical to the proof of 

(iii) Due to the length of this proof, the proof has been relegated to Appendix [B| 

(iv) The proof is identical to the proof of (iii). 


□ 
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We have established the limiting behavior of ratios of compensation parameters, coefficients, 
eigenvectors and the vector on the horizontal axis. As we have seen in Lemma 5.15, the limiting 
values of the eigenvectors i+(a,/3) and i_(a,/3) are finite. So, we can bound the absolute value 
of both of them by a constant not depending on a or (5. In doing so, one does not need to take 
into account the eigenvectors i+(a,/3) and i_(a,/3) to establish absolute convergence. Based 
on this observation and the asymptotic results derived above, we now show that the series 
appearing in (5.65) are absolutely convergent. 

Theorem 5.17. There exists a positive integer N such that 

(i) The series 


{s+iY 


s+l 


E E I“miEI’iw«)+;||/5, 

1=0 i=l j=l 

converges for all m > 0, |n| > 1 with m + |n| > N. 


H I 


(5.85) 


V+i 


(ii) The series 

oo (s+l) 

1=0 i=l 

converges for all m > 0, |n| > 1 with m + |n| > N. 


hi I 

l,i I 


(5.86) 


(hi) The series 

converges for all m > N. 
(iv) The series 


oo (s + l)* 

E E 

1=0 1=1 


(5.87) 


p{m,n,r), r = 0,1, 

m+|n|>W 


,S - 1 


(5.88) 


converges absolutely, where p{m,n,r) is given in (5.65). 

Proof, (i) We can view the series as an infinite tree with s + l roots and each term has s + l 
children. We define the ratios of a term with each of it’s descendants as 






__ |l?/+l,d(d(i)+i)+fcll«)!j)i_d(i)+jll^l+l,(i(d(i)+j)+fc 


\'ni,d(i)+j I 


a 


ohi I 


(5.89) 


We multiply the ratio j^^,{rn,n) by 


„ + l \n\ ffn+\n\ 

'^l,d{i)+j ^l + l,d{i)+i ^l,i ei,d(i)+j 
d(i) + i r.l’*! 0,1"! fl'^+l'il 


, yielding 


RlZ^ki^,n) = 


m+|n| o|/t| in 

'ni+l,d{d(i)+j)+k ^l,d{i)+j ^l+l,d(i)+j l^l,d{i)+j l^l+l,d{d{i)+j)+k ^l,i 


m+\n\ o|n| 


„hl 


,j™;h , 


a'”' B 

++i,d(i)+i Pi 


|n| 

'l,d{i)+j 


(5.90) 


5.15 


As Z —>■ oo we obtain by Lemmas 
the inequality is element-wise, and" 


and 


5.16 


that lim;^oo n) < RfY^m, n), where 


RfXk^m^n) = < 


\-in\\lv\\v-/v+\^+\'^\, j = l,2,...,s, A: = l,2,...,s, 

|6'r,||6'i..+il|ii-|'’"'k-/iii+r|l/'»^+|'”', j = s + l, A; = l,2 ,...,s, 
\- 1 r),+^\'^u\\v-/v+\^\l/v+\\'^\\w.\\'^\ j = l,2,...,s, A: = s + 1, 

j = s + l, k = s + l. 


(5.91) 
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The convergence of the series is determined by the spectral radius of the corresponding matrix 
of ratios n) := n))j^k=i, 2 ,...,s+i- If the spectral radius of the matrix n) 

is less than one, the series converges. For more details, see e.g. O Section 9] or [H proof of 

Theorem 5.2], Observe that |t;_|, |l/u+|, |rc_|, |l/t(;+| < 1, but the values of the limiting ratios 
\lv\-, |7r;s+i|, |^r?s+il b® larger than one. The spectral radius of the matrix 

R^^\m, n) depends on m + \n\ and thus we can hnd an integer N with m + n > At for which 

the spectral radius is less than one, ensuring the series converges for m + \n\ > N. 

(ii) Using a similar analysis as in (i), we define the ratio as 




Wl+2,d{d{i)+j)+k\\ci]X2,d{d{i)+j)+k\\l^\+l,d(dii)+j)+k 


10+l,d(i)+j 11 Ct]^i^ci{i)+j 11 ^1,J(i)+j I 

For 1 —)• oo we have that lim/^oo R] i 


(5.92) 


Rflki'^^^) = ^ 


i7,ii7.ii^-Kr+'”', 

16*,, 117,. 11 u_ I l”l I u_/u+1 ™ 11/t(;+11, 
l7»?.+i 11 11Ik-1I 


j = 1,2 ,... ,s, k = l,2,...,s, 
j = s + l, A; = 1,2,... ,s, 
j = 1,2,... ,s, k = s + l, 
j = s + 1, k = s + 1. 


(5.93) 


\n\ 


The spectral radius of the matrix R^‘^\m,n) := {Rj\^{m,n))j^k=i, 2 ,...,s+i depends on m + 
and thus we can find an integer N with m + n > N for which the spectral radius is less than 
one. 


(hi) We can rewrite (5.87) as 


OO (^+1) 


{s+lf 


E E i“Siii'Mi = Kiiiho,ii+E E 


1=0 i=l 


|h, 


Li I 


1=1 i=l 


+l,i\ 


(5.94) 


Lemmas 5T^iii)(a) and (iv)(a) show that —)• 0 and thus we can bound it from above 

and need not consider it when proving convergence of the series. Proving convergence of 


OO (s+l) 


! + l 


1=0 i=l 




(5.95) 


establishes convergence of (5.87). Exploiting the similarity with the series (5.86), we define 
the ratio R\ R^^{rn) := R\ Hence, the limiting ratios can be bound from above: 

lim/^oo < R^^2,.j^{m) = i?j^^(m,0). The spectral radius of the matrix R^^^m) := 

{Rj_^j^{m))j^k=i, 2 ,...,s+i depends on m and thus we can find an integer N with m> N for which 
the spectral radius is less than one. 

(iv) This follows straightforwardly, along the same lines as in [5l Section 9]. □ 

Remark 5.18. We note that the series in Theorem 5.17K i) (without the absolute values) 
corresponds to the sum of the hrst series in ( |5.65a ) and the first series in (5.65c). The series in 
Theorem 5.17[ii) (without the absolute values) corresponds to the sum of the second series in 


(5.65a) and the second series in (5.65c). So, proving convergence of the series in Theorem 5.17 
establishes convergence of the series in (5.65). 
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s 

p 

2 

0.1 

0.3 

0.5 

0.7 

0.9 

5 

0.1 

0.3 

0.5 

0.7 

0.9 

N 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


Table 2: The index N for fixed q = 0.4 and varying s and p. 


Remark 5.19. The index N is the minimal non-negative integer for which the spectral radii of 
the matrices n) and n) are both less than one for m+ |n| > N and the spectral 

radius of R^^\m) is less than one for m > N. So, for all states {m,n) with m + |n| > N and 
{N,0) the series in (5.65) converges. In general, the index N is small. In Table we list the 
index N for hxed q, while varying the values of s and p. Note that the area of convergence 
might be bigger than the one based on N, since N is based on R^^\m,n), R^‘^\m,n) and 
R^^\m) which are upper bounds on the rate of convergence. 

We are now in the position to formulate the main result of the paper. 

Theorem 5.20. For all states (m, n, r), m G Nq, n G Z, r = 0,1,..., s — 1 and m + |n| > 
N including m = N, n = 0, the equilibrium probabilities p(m, n) are proportional, up to a 


multiplicative constant C, to the expressions in (5.65), see (3.3), where C is the normalization 
constant of the equilibrium distribution. The remaining p(m, n), m + |n| <N are determined 
by the finite system of equilibrium equations for the states {m,n) with m + |n| < N, where N 
is the minimal non-negative integer for which the spectral radii of the matrices R^^\m,n) and 
R^‘^\m,n) are both less than one for m + |n| > N and the spectral radius of R^^\m) is less 
than one for m > N. 

Proof. This proof is similar to the proof in [21 proof of Theorem 5.3] and |5l Section 11], 
but we include it here for completeness. Define Cn = {(m, n) | m G No, n G Z,m + |n| > 
A^} U {{N,0,s — 1)} and note the similarity with the set of states defined in the proof of 
Lemma [5.1] and the associated transition rate diagram in Figure [l0| Then Tw is a set of states 
for which the series in (5.65) converge absolutely. The restricted stochastic process on the set 


Tat is an irreducible Markov process, whose associated equilibrium equations are identical to 
the equations of the original unrestricted process on the set Cjy, expect for the equilibrium 
equation of state {N, 0, s — 1). Hence, the process restricted to Cn is ergodic so that the series 


in (5.65) can be normalized to produce the equilibrium distribution of the restricted process 
on Cn- Since the set Nq x Z x {0,1,..., s — 1} \ Cn is hnite, it follows that the original 
process is ergodic and relating appropriately the equilibrium probabilities of the unrestricted 
and restricted process completes the proof. □ 

The following remark is in line with Remark |4.1[ 

Remark 5.21. {An efficient numerical scheme) The following scheme exploits the fact that 


the rate of convergence of the series in (5.65) increases as m + |n| increases. So, further away 


from the origin, fewer compensation steps L are needed to achieve the same accuracy according 
to (4.1) and (4.2). This property was seen in Section]^ and in particular Figure Denote a 


triangular set of states as Tx '■= {{m,n) 
visual aid. 


m G No, n G Z, m + |n| < x}. Figure 12 


serves as a 


(i) Determine the minimal non-negative integer N for which the spectral radii of the matrices 
R^^^ (m, n) and (m, n) are both less than one for m -|- |n| > N and the spectral radius 
of R^^'^ (m) is less than one for m > N. 

(ii) Select integers M and K such that N < M < K. 
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(iii) Determine p(m, n), (m, n) G Tk \ Tm according to (4.1) with (7 = 1 and L the minimal 
integer such that (4.2) holds. 


(iv) Determine p(m,n), {m,n) G Tm from the equilibrium equations for the states {m,n) G 
Tm- 


(v) Normalize the equilibrium distribution by dividing each equilibrium probability by the 
sum E(m,n)erKP("^’”)l- 

The integer L in step (iii) depends on M. As M increases, L decreases or stays constant, but 
the size of the system of equilibrium equations in step (iv) increases. This tradeoff is clearly 
in favor of selecting a larger M: decreasing L decreases the number of computation steps 
exponentially, whereas the size of the system of equilibrium equations increases polynomially 
with M. As K increases, the equilibrium probabilities become more accurate. K can be chosen 
arbitrarily large; it has little to no impact on the performance. 


6 Conclusion 

We have studied a queueing system with two non-identical servers with service rates 1 and 
s G N, respectively, Poisson arrivals and the SED routing policy. This policy assigns an 
arriving customer to the queue with the smallest expected delay, i.e. waiting time plus service 
time. The SED routing policy is a natural and simple routing policy that balances the load for 
the two non-identical servers. Although not always optimal, SED performs well at both ends of 
system utilization range. Moreover, SED routing is asymptotically optimal in the heavy traffic 
regime. 
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The SED system can be modeled as an inhomogeneous random walk in the quadrant. By 
appropriately transforming the state space, we mapped the two-dimensional state space into 
a half-plane with a finite third dimension. The random walks on each quadrant are different, 
yet homogeneous inside each quadrant. Extending the compensation approach to this three- 
dimensional setting, we showed in this paper that the equilibrium distribution of the joint queue 
length can be represented as a series of product-form solutions. These product-form solutions 
are generated iteratively to compensate for the error introduced by its preceding product-form 
term. 

The analysis presented in this paper proves that the compensation approach can be applied 
in the context of a three-dimensional state space. We believe that a similar analysis can be 
used to investigate general conditions for applicability of the compensation approach for three- 
dimensional Markov processes. These conditions will be comparable to the three conditions in 
Section[^ Furthermore, the compensation approach can possibly be extended to the case where 
all three dimensions are infinite, paving the way for performance analysis of higher-dimensional 
Markov processes. 

The insights gained for the SED system with two servers, specifically, the series expressions 
for the equilibrium probabilities, can be used to develop approximations of the performance of 
heterogeneous multi-server systems with a SED routing protocol. These approximations can 
be derived along the same lines as in [2T]. An approximate performance analysis for a system 
with two servers can be found in |28| . 

Another interesting direction for future research is to study rare events or tail probabilities 
in the SED system, in a similar way as done for JSQ systems in m Since the compen¬ 
sation approach determines the complete expansion of each equilibrium probability, one can 
approximate rare events with arbitrary precision. 
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A Proof of step (c) of Lemma 5.6Ki) 


The following argument was communicated to us by A.J.E.M. Janssen. Our goal is to show 
that for every a with |a| G (0,1) the equation 


al3{l + s){p -I-1) - /3^(1 -I- s)p - 
a 13 s 


= i = l,2 ,...,, 


(A.l) 


has at least one root Pi with |/Jj| < |a| for each i. To that end, we consider the more general 
formulation 

a = f{z) = ^{v+-z){z-v_) (A.2) 

with z = P/a, E = > Q, t = ^ G (1, 2] and v± is defined in (5.27) with the properties 

1 


0 < < 1, u+> 1, v+ + v_ = ^^^, v+v_ = ———. 

p (l + s)/3 


(A.3) 


One retrieves the original form (A.l) by setting a = Uio/^!^. 
The proof consists of three main steps: 
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(a) We derive the inverse function of a = f{z) on a neighborhood of z = v_ using the 
Lagrange inversion theorem m Section 2.2], We establish that z = g{a) where g{a) is a 
power series that is analytic on a neighbourhood of u = 0 and has positive power series 
coefficients. 


(b) Employing Pringsheim’s theorem [8l Theorem 1], we are able to extend the radius of 
convergence for the power series g{a) to Umax > 1- 

(c) Finally, we establish that for a = i = 1,2,..., s satisfying (A.2), the correspond¬ 

ing 2 ; has \z\ < 1. 


(a) We note that f{z) is analytic near z = v_. So, by the Lagrange inversion theorem, we have 
in a neighborhood of u = f{v_) = 0 


2 ^ = 5(<7) =v_ + 


V —— 

^ n\ dz”“^ ^ 

n=l 
00 


Z — V_ 


= V, 


+ E 


^iv+- z){z -v_) 

{a/EY d”"^ 


n=l 


n! Az'^~^ — zY 


(A.4) 


Note that g{-) is analytic in a neighborhood of u = f{v_) = 0. Let gi{z) = z"’* and g 2 {z) = 
(u+ — 2 :)“"’ so that we can write 


4 n —1 


n—1 


dz' 


i::i9i{z)g2{z) = Yj 


k=0 


n — 1 


^n-k-l ^k 


where 


^n — h—^ n—k—2 

" =(n 


dz 


-z z 


nt—(n—k—l) 


2=0 


— n^(z) = (^ + ^~ 0 ' ^, _ ^)-(n+k) 


(A.5) 


(A.6) 

(A.7) 


Thus, we can conclude that 


d 


n—1 


^nt 


d2;n-l 


n-1 , , ..N, n-k-2 

^-■v- k=0 ^ ' i=0 


— I 


nt—{n—k—l) 




/ (v. 

_ 'fj ^ 21-|- k 


> 0 . 


(A.8) 


Hence, the power series g{cr) has positive power series coefficients. 

(b) We now extend the convergence range of g{a) from a neighbourhood of 0 to the entire 
unit disc. This is achieved by using the Pringsheim theorem. We are allowed to do so because 
the coefficients of the power series of g{cr) are positive. According to the Pringsheim theorem 


we can limit attention to the positive real line. With reference to Figure 13, we let Umax = 
max{/( 2 ;) | 2 ; > u_}. Note that for 2 ; G [u.,^(umax)) the function f{z) is strictly increasing and 
is analytic by ( A.2| ), thus its inverse, g{a), is analytic for u G [/(i'-),/(fl'(o'max))) = [0. U max ). 
Clearly, Umax > Ij hence by the Pringsheim theorem g{a) is analytic for |u| < Umax and thus 
inside the unit disc. 
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(c) It is important to see where g{(7raax) lies. We find z = g(a max ) from 


0 = 


(v+v_t -{t- l)(r;+ + v_)z - (2 - 

{TT^ - (^ - 1 )—^ - (2 - t)zA . 
z*+iV(l + s)p ^ ’ p ^ W 


(A.9) 


For t G (1, 2) we take the positive solution of the quadratic equation in (A.9) and obtain after 
some rewriting 


^ — fl'(crmax) — 


1 + p + ^/{l + pY + 4p(s - 1) 


< 1 , 


(A.IO) 


which also holds for t = 2. So, for a cr with |it| < cimax, we have from the positivity of the 
power series coefficients that 


|5'(cr)l < g{W\) < 9(o-max) < 1- 


(A.ll) 


By taking cjj = z = 1, 2,..., s with \ai\ < 1 we see that the power series g{(Ti) converges 

and l^l = \g{ai)\ < 1. Finally, we can conclude that for every a with |a| G (0,1), (ig has at 
least one root Pi with |/3j| < |q;| for each i. 


B Proof of Lemma 5.16(iii) 


In this appendix we prove the convergence of the ratio of coefficients used in the horizontal 
compensation step for a solution that satisfies the positive inner equations. Specifically, we 
consider the setting of Lemma 5.12[ i). We wish to establish the limiting values, or an upper 
bound on the absolute value of the following ratios, as a 0, 


(a) h/i/, 

(b) Ps+i/y, 

(c) pi/v, z = 1,2 ,... ,s. 
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We reiterate below the equations that determine these ratios, which can be found in the main 
body of the paper in (5.58)-(5.60), 


(^0,1 + a^-i,i)h - a ^ r]d+{a, f3i) = iyai+{a, /3), 


(B.l) 


i=l 


aBQflh + _i + aAo _i)i+(a,/3i) 


2=1 


+7ys+i/3s+i(-Bi4 + aBo,i)i-{a, /3s+i) = + a^o-i)i+(a,/3), (B.2) 

—rjs+ias + ash{0) + (1 + s)pqh{s — 1) = 0. (B.3) 


As it turns out, when we let a ^ 0 in the system of linear equations (B.1)-(B.3), we obtain a 
degenerate system of equations. By properly scaling the system, one obtains a non-degenerate 
system of equations. 

In this appendix we adopt the following notation to indicate the limiting value of a ratio 


:= lim - 
QtO V 


X 


(B.4) 


(a) Divide (B.3) by u and rearrange to obtain 


h{s-l) s 1 fns+i h{0) 
-= a- 


. , , , ,. (B.5) 

u 9(1 + s)p \ ly V J 

By letting a J, 0 in (B.5) we get — 1) = 0. Next, divide (B.2) by a and v, which yields 

Bofl —h ^ — — (Ai 1 + Q!Ao,-i)i+(a, A) 

7/ z—/ ' 


V V OL 

2 = 1 

+ + uBq I3s+i) = --{M-i + aAo i)\+{a, 13). 

u a a 


(B.6) 


We let a i 0 and use that Ai^-i = (1 + s)pl and Bip = (1 + s)pM^^'^ This gives 

S 

Bo,oh^™ + u_(l + s)p^r/|™eo + + s)p'ip+iO)^~^eo = -u+(l + s)peo. (B.7) 


2=1 


The matrix Bq^ is a tri-diagonal matrix and together with — 1) = 0 we find that = 0. 

(b) Together with = 0, equation (B.7) reduces to the single equation 




(B.8) 


2=1 


Now, divide (B.l) by a and u, yielding 


h h * 

Ao,i-h A_i4- y^ nii+{a, Pi) = i+(a,/3). 

au u 

2=1 


(B.9) 


Note that Ao,ih/(ai/) = (1 + s)/9(l — q)h{^s — l)/(ai/)eo, so that, together with (B.5) and 

(B.IO) 


A_i^i = /, (B.9) reduces to 

1-9/^s+i ^( 0 ) 


q \ V 


) Xl. 

eo + - - A) = i+(a,/3)- 


2 = 1 
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Moreover, letting a | 0 again gives us a single equation, namely 




(B.ll) 


2=1 


One can solve the linear system of equations (B.8) and (B.ll) for the two unknowns 
and where the solution of is given in Lemma |5.16 'iii)(b). 

(c) We wish to establish a scaling of each equation in ( B.1[ )-(B.3 ) so that, in the limit for 
a 0, we obtain a non-degenerate system of equations. Fr om (5.14) and (5.25) we know that 
i+{a, I3i, r) = ,r = 0,1,..., s — 1 and from Remark 5.7 it is clear that there exists a j 

for which z+(a,/3,r) = Uj(3''^^, r = 0,l,...,s — 1. Thus, the correct scaling to establish a 
non-degenerate limit of a term i_,_(a,/3,r) is 

We are going to multiply both ( |B.6 ) and (B.9) element-wise by the column vector 

Q := (l a.-^/s a-2/s ... 

To be able to do that, we need some new notation. Let us introduce the column vectors 


. _ (m Mil 


h{s—l) 


V ■= {rji rj 2 ■■■ Vs) , and h :: 
and the matrices 

W(a) := (i+(a,/3i) i+(a,/32) ••• i+(a,^s)), 

^i+(a,/32) ••• /?.))■ 

Using the introduced vectors and matrices we can write 


'^—i+{a,/3i) = W{a)-, and ^ —i+(a,/3j) = W(a)-. 

V V V a V 

2=1 2 = 1 


(B.13) 

(B.14) 

(B.15) 

(B.16) 


We introduce the symbol o for element-wise (Hadamard) multiplication, i.e., for column vectors 
X and y, we have that z = xoy means that element z{r) = x{r)y{r). Below we list some useful 
element-wise multiplications that we will use 




/ 1 


1 


1 \ 



W{a) = 


^02 




, (B.17) 

a o 









(&)(. 





-.OLO 

(a,/3) = 1 


■■ 0“ 


T 


(B.18) 


Q o h = 





(B.19) 

Moreover, for 

a 1 0, we 

determine that 

OLoW (a) 

— )• W, OL O 

W{a 

—)> V-W, and a 

oi+{a,/3) 


Wj, where W and Wj are given in (5.78)-(5.79). 

We are now in a position to multiply both (B.6) and (B.9) element-wise by a. For the 
element-wise multiplication of (B.6) we get 

« o ° E -- (^1,-1 + A) + o + aBo,i)i_(a, /3.+i) 

7/ ^ ^ V(y.^ ' 1/ n/ ^ ' 


2 = 1 


V a 




= --a o (^ 1-1 a^o-i)i+(a, /3)- 


a 


(B.20) 
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Combining (B.20) with the definition of the transition matrices i?o,o = (1 + s)pL — (1 + s){p + 
1)1 + sL?", ^ 1,-1 = (1 + s)pl and Bi^i = (1 + s)pM^^’^~^\ and the matrix-vector notation 
described above, we obtain 

^(1 -|- s)pL — -|- s){p + 1)1 + (1 -|- s)poL o Ty(a) — 


-|- Q; o 


E 

2=1 


Vil^i . ?/s+l l^s+l 

- aA{) Pi) H- 

V a V a 


= --OL O ((1 -h s)pl + aAo i)\+{a,l5). 
a ' 

Multiplying (|B.9|) element-wise by cx gives 


((1 -F s)pM^^'^ -h ao: o Bo,i)L(a,/3s+i) 

(B.21) 


^(1-|-s)p(l — — — a o Ty(Q;)—= a o /3). (B.22) 


Finally, if we let a ^ 0 in (B.21) and (B.22) we obtain the limiting system of linear equations 

Lh'™ -h = -v+Wj - 7^^_^it(;_'(/;+(0)^"^eo, (B.23) 

(1 + s)p(l - - -Fry'™ = w^. (B.24) 

We have constructed 2s equations, (B.23)-(B.24), for the 2s unknowns and 77 *™. The 
solution to this system of equations depends on j. So, let h*™ and rjp be the solution to 
(B.23)-(B.24) for a specific j. Then, 


lim 

atO 


<^^max^ \pP{r)\. 


l^j,r+l<s 


(B.25) 


This concludes the proof of part (c). 
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